The Effect of Systematic Error in Forced Oscillation Testing 
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One of the fundamental problems in flight dynamics is the formulation of aerodynamic forces and 
moments acting on an aircraft in arbitrary motion. Classically, conventional stability derivatives are used for 
the representation of aerodynamic loads in the aircraft equations of motion. However, for modern aircraft 
with highly nonlinear and unsteady aerodynamic characteristics undergoing maneuvers at high angle of 
attack and/or angular rates the conventional stability derivative model is no longer valid. Attempts to 
formulate aerodynamic model equations with unsteady terms are based on several different wind tunnel 
techniques: for example, captive, wind tunnel single degree -of-freedom, and wind tunnel free-flying 
techniques. One of the most common techniques is forced oscillation testing. However, the forced oscillation 
testing method does not address the systematic and systematic correlation errors from the test apparatus that 
cause inconsistencies in the measured oscillatory stability derivatives. The primary objective of this study is 
to identify the possible sources and magnitude of systematic error in representative dynamic test 
apparatuses. Sensitivities of the longitudinal stability derivatives to systematic errors are computed, using a 
high fidelity simulation of a forced oscillation test rig, and assessed using both Design of Experiments and 
Monte Carlo methods. 
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Nomenclature 

damping constant [N*m/(rad/s)] 
backlash [arcmin] 
mean aerodynamic chord [m] 
lift coefficient 

lift coefficient due to angle of attack 

change in lift coefficient due to angle-of-attack rate 

change in lift coefficient due to pitch rate 

change in lift coefficient due to pitch acceleration 

change in pitching moment coefficient due to angle of attack 

change in pitching moment coefficient due to angle-of-attack rate 

change in pitching moment coefficient due to pitch rate 
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change in pitching moment coefficient due to pitch acceleration 

diagonal elements of "hat” matrix (see ref. 32) 

input saturation 

moment of inertia [kg*m 2 ] 

reference length 

reduced frequency, k=col/V 

predicted error sum of squares 

pitch angular velocity [rad/sec] 

coefficient of determination 

sum of squares 

error sum of squares 

total sum of squares 

regression sum of squares 

time [sec] 

period [sec] 

motor torque [N*m] 

load torque [N*m] 

velocity [m/s] 

increment 

angle of attack [rad] 

oscillation angle of attack (amplitude) 

regression coefficients 

increment 

error 

pitch angle [rad] 
angular velocity [rad/sec] 


Superscripts 

= oscillatory data (in-phase or out-of-phase) 

Subscripts 

adj = adjusted 

EQ = equivalent 


I. Introduction and Problem Definition 

W ITH the increasing need for high maneuverability capabilities in modern combat aircraft and for improved 
safety of transport aircraft under loss of control conditions there are significant shortcomings in using the 
conventional stability derivative modeling approach. Current prediction methods work adequately over low to 
moderate angles of attack. However, at high angles of attack (pre- and post-stall conditions) nonlinear unsteady 
aerodynamics are more pronounced. In this region the aerodynamic characteristics are based on shock waves, 
separated flows, vortical flows, movement of separation points, and variation of vortex breakdown locations over 
the surface. 1 

Currently, there are two primary approaches taken to formulate the nonlinear unsteady aerodynamic models - (1) 
analytical techniques and (2) experimental techniques. Current analytical approaches to modeling the unsteady 
aerodynamic characteristics of an aircraft vary greatly, and there appears to be no standard methodology. The 
various analytical techniques used have been reviewed by Kyle et al. and Greenwell. 12 Dynamic test techniques 
range from captive, wind tunnel single degree-of-freedom and free-flying, and outside free-flying. Additional wind 
tunnel dynamic test techniques have been reviewed by Owens et al. and Tomek et al. 3 ' 4 The present study focuses 
on the most commonly used method, forced oscillation testing in a wind tunnel. Forced oscillation experiments are 
typically conducted at small amplitudes (< 5 deg) and single frequencies ranging from 0.1 to 10 Hz. An example of 
a forced oscillation system (FOS) in the NASA Langley 12-ft Low Speed Wind Tunnel is shown in Figure 1. 
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Figure 1. FOS installed in the 12-ft Wind Tunnel (roll configuration) 


There are several limitations in forced oscillation testing. The apparatus used for forced oscillation is often large 
in size compared to the model and in close proximity to the model, causing support inference concerns. Other 
concerns lie in the accuracy of the force and moment measurements caused by bias and random errors from 
instrumentation, calibrations and tare interpolations, structural dynamics, vortex shedding, and wind tunnel fan beat 
frequencies. *’ 514 The studies of Kim and Heim clearly demonstrated the need for proper calibration techniques and 
the need for adequate data sampling; both studies assumed that systematic errors could be removed by calibration 
alone . 15 16 Neither study looked at the dynamic effects of the test apparatus interfering with the aerodynamic 
measurements. 

Understanding the dynamics of a test rig is important because of the possible influence of the nonlinear behavior 
of components of the test rig on an aerodynamic measurement. The interactions between the mechanical 
components, transmission components, and electrical components introduce systematic errors that can affect the 
aerodynamic measurement. If the level of nonlinearities in a system can be neglected without exceeding the error 
tolerance then the system can be assumed linear. However, this assumption is not valid if the system nonlinearities 
are seen to interact with the aerodynamic loads, which pose a problem for accurately obtaining the dynamic stability 
derivatives that are used to develop aircraft control laws and flight dynamics studies. There is significant concern as 
to the overall accuracy of the measurement, particularly with forced oscillation experiments conducted at very low 
frequency ranges. This study attempts to examine the level of the systematic error in dynamic stability derivatives 
caused by mechanical and transmission components in forced oscillation test apparatuses using a high fidelity 
computer simulation of a dynamic pitch test rig, and assessed using the fundamental principles from Design of 
Experiments (DOE), and Monte Carlo methods. The study is limited to longitudinal motion only and assumes that 
mechanical vibrations are minimal and that aeroelastic effects may be negated. 


II. Theory 


A. Determination of Stability Derivatives from Experimental Data 

Dynamic longitudinal stability derivatives cannot be measured directly because the angle-of-attack rate and pitch 
rate are physically identical during a forced oscillation test and it is not possible to separate the components 
associated with pitch rate and a. Therefore, the stability derivatives will take the form of a combination; for 
example C m& + C mq . The aerodynamic coefficients are assumed to be linear functions of the angle of attack, pitching 
velocity, and their associated rates. The incremental lift coefficient with respect to its mean value is formulated as: 


A C L = C L Aa + 


C , cc + 

‘-‘n 


C L q + 


O' 

vVO 


C,j 


(1) 


For harmonic motion. 
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( 2 ) 


Then, 


A a = a A sin cat 
a = q = coa A cos cat 
a = q = -ca 2 a A sin cat 

A C L = a A ( C L -k 2 C L )sin cat + a A k(c L + C L )cos cat 
= a A (c, sin cat + kC L cos cat ) 


(3) 


where the in-phase and out-of-phase components of lift coefficient have been identified as C L and C L , 

respectively. Using the same approach for pitching moment, the in-phase and out-of-phase pitching moment 
coefficients are, respectively. 


c 

c 


= C„ 

= c 



(4) 


The in-phase and out-of-phase components of lift coefficient are found by integrating the time histories of ACl over 
a given number of cycles as: 


n r T 


C L = 


C L = 


a A n c T 0 


}AC t (r) 


sin cat dt 


a A n c kT 0 


jAC t (t) 


cos cat dt 


(5) 


where T = 2n/(» is the period. The same approach is used to compute the in-phase and out-of-phase components for 
pitching moment, so that: 

sin cot dt 


C = — f AC m (f)cos cot dt 

q a A n c kT J 0 

Equations (5) and (6) have assumed that the sine waveform is of high fidelity. 


nj 


c = 


a A n c T 0 




B. NASA Langley 12-ft Wind Tunnel Dynamic Test Rig 

A new forced oscillation system (FOS) was developed at the NASA Langley 12-ft Low Speed Wind Tunnel to 
improve operational efficiencies. The wind tunnel can operate at dynamic pressures up to 7 psf ( U , = 77 ft/sec at 
sea-level conditions) which corresponds to a unit Reynolds number of approximately 492,000 per foot. The test 
section has a turbulence level of about 0.6 percent. 17 The test section airflow is produced by a 15.8-ft diameter, 6 
blade drive fan powered by a 280 HP, 600 V, 600 rpm DC motor and controlled by a 500 Hp AC motor. 17 

The system supports both static and dynamic testing, and is capable of independent setting of a and (1 
combinations. For pitch oscillation testing, the system has an angle-of-attack range from -5° to 85°. The maximum 
pitch capability of the test rig is 260 °/sec pitch rate and 2290 °/sec 2 pitch acceleration. For roll oscillation testing, 
the system has a range of displacements from ±170 degrees with a maximum capability of 190 °/sec roll rate and 
12,750°/sec 2 roll acceleration. The pitch and roll capabilities of the test rig are the no-load limits. Figure 2 
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illustrates the overall FOS mechanical subsystem. The user prescribes the position motion, and the control 
algorithm executes a positional, closed-loop algorithm for tracking. The velocity control loop is managed internally 
using a power (BDS4) amplifier. The amplifier powers the motor which imparts the desired test article motion via 
the drive system components. 


Electronics and Control Hardware FOS Tunnel Hardware 

(Located outside test section) (Located in test section) 



Figure 2. Overall FOS block diagram. 


C. Drive Motor and Compliantly-Coupled Drivetrain Model 

The Kollmorgen motor is a three-phase AC induction motor. A detailed derivation and explanation of the 
equations of motion for the motor can be found in reference 30. The desire for high performance motor drives is 
necessary in applications that require rapid command response and dynamic stiffness. The downside of high 
performance is mechanical resonance which is caused by compliance between two or more components in the 
mechanical transmission chain. For the case being studied, the resonance is typically the compliance between the 
motor and load. 

Resonance can occur in a forced oscillation test when operating near resonant frequencies of the drivetrain. For 
the system considered in this study instability was unlikely with test frequencies between 0.1 and 10 Hz and the 
spring constants very large. Further analysis proved that the compliantly-coupled drivetrain can be modeled as a 
rigid body and the motor torque equation becomes simply. 


T 

motor 




( B bal + B sting) + 




(7) 


The T load represents the aerodynamic load from the test article and the gear reduction ratio is given as a ratio of the 
number of teeth (N). 


III. Computer Simulation Methodology 

A computer simulation of the FOS was created using a simplified model. Aerodynamic loads for the model 
were obtained by using regression models built from an existing experimental data set. The simulation was then 
used to perform a sensitivity study subject to changes in sting inertia, system damping, backlash, reduced frequency 
and input saturation. 

A. Computer Simulation Details 

Figure 3 provides a conceptual block diagram of the computer simulation created in Simulink©. Starting from 
the left, the user commands a forced oscillation profile (typically a sine wave). The commanded motion is 
controlled through a proportional-integral-derivative (PID) controller and a velocity signal is output and filtered 
before the command signal is amplified. The velocity signal is then fed into another PID controller, filtered, and 
saturated. Normally, a microprocessor monitors the incoming signal and then controls the motor. This step could 
not be modeled in Simulink due to the unknown internal process of the microprocessor. However, in this simulation 
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after the saturation of the velocity signal, a three-phase current signal is modeled and output. The incoming current 
signal is controlled by a PI controller. The output is then transferred to a pulse-width modulator which controls and 
powers the motor. The torque of the motor is output, travels through the mechanical drive train and compliantly 
coupled system to actuate the aerodynamic model. Current, velocity, and position are fed back to their respective 
PID controllers 30 . 


Plant 



To represent the aerodynamic loads, the study implemented an experimental model data set due largely to the 
lack of viable Computational Fluid Dynamics (CFD) codes for predicting high angle of attack and unsteady 
aerodynamic effects for an entire aircraft. The data set was taken from NASA LaRC 12-ft Low Speed Wind Tunnel 
experiments during the 1990s. 18 The chosen test article was a 10% scale model of an F16-XL aircraft. An ordinary 
least squares regression model was fit to the data as shown below to define the in-phase and out-of-phase lift and 
pitching moment coefficients in the simulation. The in-phase and out-of-phase coefficients are a function of both 
angle of attack and reduced frequency. Regression models were developed at a fixed angle of attack with reduced 
frequency as a variable. As an example, equations 8-11 show the models built for an angle of attack of 30.8°. 


C, = 82.2k 3 -15.4k 2 +23. 3k- 0.387 (8) 

L -‘a 

C, = -l\28k 3 +1005& 2 -298& + 31.6 (9) 

C m = 6.82k 3 - 8 Aik 2 +3.19/: -0.0191 (10) 

rn a 

C m = -21 4k 3 + 228k 2 -62.9k + 5.48 (11) 

m q 


The dynamic contribution was determined by Eq. (12) and (13) after subtracting the static contributions. 


m dyn 


= c L o+ 

c 

Cj 0 

(12) 


v2 V, 

L q 


( r 3 



= c m e+ 

v 2V y 

C ’+° 

(13) 


Where the forcing function is based on harmonic motion, see Eq. (14). 
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(14) 


6 ~a A sin cot 
0 = a A co cos cot 

The stability derivatives are computed using numerical integration of equations 5-6 and stored to a file. 

B. Computer Simulation Model Verification and Validation 

The computer simulation was validated and verified by using several steps. The first step was to insure that the 
theories and assumptions underlying the conceptual model are correct. The second step was to certify that the 
computerized model reflects the physical system based on the conceptual model. The final step was to validate 
individual subcomponents (i.e. motor model and aerodynamic model) of the system against either historical data or 
benchmark cases. An example is provided in Figure 4. Here the in-phase lift coefficient predictions are compared to 
the experimental data over a range of reduced frequencies and angles of attack. It should be noted that although the 
submodel was validated against various reduced frequencies over a range of angle of attack, only the regression 
models for a fixed angle of attack of 30.8° and variable reduced frequency (refer to Eq. 8-11) were used in the final 
simulation. These regression models were chosen as representative of the general characteristics associated with the 
unsteady aerodynamics. For further details the reader should refer to Williams. 30 



20 25 30 35 40 45 50 55 60 65 

a [deg] 


Figure 4. In-phase lift coefficient validation. 

IV. Experiment Design 

The design chosen for this study is called a nested central composite face centered design and is based on the use 
of two, 2 V 5 "’ fractional factorials augmented with face centered axial and center points. 19 The design points of the 
Nested FCD in two-factor space are illustrated graphically in Figure 5. The extremes of the factor levels are set on 
the perimeter of the outer box. The nested factor levels are set on the perimeter of the inner box. The design 
supports regression modeling with pure cubic terms in addition to the full quadratic model. It should be stated that 
face centered designs are not rotatable. However, it is generally not a priority when the region of interest is 
cubeoidal. 20 Adding one or two center points is sufficient to produce reasonable stability in the variance of 
prediction. 
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outer factorial point 


Figure 5. Nested face centered design in two factors. 19 


Table 1. Factor limits. 


Factor 

Units 

Outer Low 

Inner Low 

Center 

Inner High 

Outer High 

A 

kg*m 2 

0.0001729 

0.0002944 

0.0004159 

0.0005374 

0.0006589 

B 

N-m-s/rad 

0.0000000 

0.0020140 

0.0040280 

0.0060420 

0.0080560 

C 

- 

0.0810000 

0.0857500 

0.0905000 

0.0952500 

0.1000000 

D 

arcmin 

2.0000000 

4.0000000 

6.0000000 

8.0000000 

10.0000000 

E 

rad 

0.0809000 

0.0856775 

0.0840850 

0.0856775 

0.0872700 

Factor Labeling Key: 


Label 




Factor 




A 



Jeq (equivalent inertia) 




B 



Beq (equivalent damping) 




C 



k (reduced frequency) 




D 



BL (backlash) 




E 



IS (input saturation) 



The five factors investigated in this study were: equivalent inertia, equivalent damping, reduced frequency, 
backlash, and input saturation. The responses were: in-phase lift coefficient, out-of-phase lift coefficient, in-phase 
pitching moment coefficient, and out-of-phase pitching moment coefficient. The factor limits for the design are 
specified in Table 1. Equivalent inertia is based on the compliance of the system, assuming a straight sting (Outer 
Low) or a bent sting (Outer High). Equivalent damping is similar; however, it assumes no damping (Outer Low) 
and the system compliance damping (Outer High). The reduced frequency range is based on experimental data. 
The backlash range in the gear drive is based on published data from the manufacturer. 21 Finally, the input 
saturation range is based on observations from experimental data noted by Kim et al. 22 Experimental results have 
shown evidence of clipping in the data. For example, a sinusoidal wave form is input into the system and the output 
should be a sinusoidal waveform; however, experimental results showed the peaks of the waveform being clipped. 
This clipping was traced back to input saturation. The inner face center design points were determined by Vi of the 
distance to the outer face center design points. 

The nested face centered design test matrix was analyzed using Design Expert©, commercial software; see 
Table 2. The simulation was run for each test case listed in Table 2. Care must be taken when analyzing results of a 
response surface methodology with a deterministic computer simulation. 


8 


American Institute of Aeronautics and Astronautics 





Table 2. 2 5 ' 1 nested face centered design test matrix (center run is shown as bold). 


Run 

Jeq (A) 

Beq(B) 

k(C) 

BL (D) 

IS (E) 

Run 

Jeq 

Beq 

k 

BL 

IS 

1 

0.0006589 

0 

0.081 

10 

0.08727 

30 

0.0002944 

0.006042 

0.08575 

4 

0.0824925 

2 

0.0004159 

0.004028 

0.0905 

6 

0.084085 

31 

0.0005374 

0.006042 

0.08575 

4 

0.0856775 

3 

0.0004159 

0.004028 

0.0905 

10 

0.084085 

32 

0.0002944 

0.002014 

0.09525 

4 

0.0824925 

4 

0.0006589 

0.008056 

0.1 

2 

0.0809 

33 

0.0005374 

0.002014 

0.09525 

4 

0.0856775 

5 

0.0004159 

0.004028 

0.1 

6 

0.084085 

34 

0.0002944 

0.006042 

0.09525 

4 

0.0856775 

6 

0.0004159 

0.004028 

0.0905 

2 

0.084085 

35 

0.0005374 

0.006042 

0.09525 

4 

0.0824925 

7 

0.0006589 

0.008056 

0.081 

10 

0.0809 

36 

0.0002944 

0.002014 

0.08575 

8 

0.0824925 

8 

0.0001729 

0 

0.1 

10 

0.08727 

37 

0.0005374 

0.002014 

0.08575 

8 

0.0856775 

9 

0.0006589 

0 

0.081 

2 

0.0809 

38 

0.0002944 

0.006042 

0.08575 

8 

0.0856775 

10 

0.0006589 

0.004028 

0.0905 

6 

0.084085 

39 

0.0005374 

0.006042 

0.08575 

8 

0.0824925 

11 

0.0001729 

0.008056 

0.1 

10 

0.0809 

40 

0.0002944 

0.002014 

0.09525 

8 

0.0856775 

12 

0.0004159 

0 

0.0905 

6 

0.084085 

41 

0.0005374 

0.002014 

0.09525 

8 

0.0824925 

13 

0.0004159 

0.004028 

0.081 

6 

0.084085 

42 

0.0002944 

0.006042 

0.09525 

8 

0.0824925 

14 

0.0006589 

0.008056 

0.1 

10 

0.08727 

43 

0.0005374 

0.006042 

0.09525 

8 

0.0856775 

15 

0.0006589 

0 

0.1 

10 

0.0809 

44 

0.0002944 

0.004028 

0.0905 

6 

0.084085 

16 

0.0001729 

0 

0.081 

2 

0.08727 

45 

0.0005374 

0.004028 

0.0905 

6 

0.084085 

17 

0.0001729 

0.008056 

0.081 

10 

0.08727 

46 

0.0004159 

0.002014 

0.0905 

6 

0.084085 

18 

0.0001729 

0.004028 

0.0905 

6 

0.084085 

47 

0.0004159 

0.006042 

0.0905 

6 

0.084085 

19 

0.0001729 

0 

0.081 

10 

0.0809 

48 

0.0004159 

0.004028 

0.08575 

6 

0.084085 

20 

0.0001729 

0 

0.1 

2 

0.0809 

49 

0.0004159 

0.004028 

0.09525 

6 

0.084085 

21 

0.0004159 

0.004028 

0.0905 

6 

0.08727 

50 

0.0004159 

0.004028 

0.0905 

4 

0.084085 

22 

0.0006589 

0 

0.1 

2 

0.08727 

51 

0.0004159 

0.004028 

0.0905 

8 

0.084085 

23 

0.0001729 

0.008056 

0.081 

2 

0.0809 

52 

0.0004159 

0.004028 

0.0905 

6 

0.0824925 

24 

0.0004159 

0.004028 

0.0905 

6 

0.0809 

53 

0.0004159 

0.004028 

0.0905 

6 

0.0856775 

25 

0.0006589 

0.008056 

0.081 

2 

0.08727 







26 

0.0001729 

0.008056 

0.1 

2 

0.08727 







27 

0.0004159 

0.008056 

0.0905 

6 

0.084085 







28 

0.0002944 

0.002014 

0.08575 

4 

0.0856775 







29 

0.0005374 

0.002014 

0.08575 

4 

0.0824925 








A. Statistics and Deterministic Computer Models 

Response surface methodology (RSM) has been used successfully with computer simulation models of physical 
systems. A few examples are Barton (1992, 1994) and Simpson et al. (1997). 23 25 The goal of RSM is to build a 
regression model of the system being modeled by the computer simulation - a regression metamodel. Although 
RSM is primarily intended for applications with random error, deterministic applications are appropriate and the 
analysis is simplified. Care must be taken when using statistics to develop the regression metamodel. The nature of 
the computer simulation is represented as an input-output relation, 

y = f W as) 

Then an estimated output is represented as 

y = g{x)+£ hias (16) 

where g(x) is the metamodel, and s bias is the error of approximation. The relationship conflicts with the use of 
Analysis of Variance (ANOVA) with least squares regression. 26 Consequently, the response surface fitting is 
hampered by the following constraints. 27, 28 

1 . The response surface model adequacy is determined solely by systematic bias. 

2. The usual measures of uncertainty derived from least squares residuals have no obvious statistical meaning; 
deterministic measures of uncertainty exist. 
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3. The notions of experimental blocking, replication, and randomization are irrelevant. 

Statistical measures such as pure error, F-statistics, and mean squared error for verification of model adequacy have 
no statistical meaning because they are all based on random error. Care should be taken when using commercial 
software because manufacturers employ these types of statistical measures for developing a model. These measures 
have no statistical meaning since they assume that observations include an error term, which has a mean of zero and 
a non-zero standard deviation. The relative size of the coded (centered and scaled) regression model terms and sum 
of squares (SS) contributions may be compared. 


B. Fitting and Validating of Regression Metamodels 

The best metrics to verify and validate the metamodels (i.e. regression models of the in-phase and out-of-phase 
coefficients determined by the computer simulation) are the coefficients of determination, cross-validation and point 
prediction methods , 26 ' 27 The classic coefficient of determination measure is 


_ SS R 


SS T 


_ 1 _ S$E 


SS T 


where 0 < R <1 


(17) 


r J' rp 

Values of R~ that are close to one imply that the variability in the response is explained by the regression model. 29 

Caution should be exercised when using R~ . Increasing the number of regression model terms always increases R~ 
regardless of the value of the contribution of the term. This does not necessarily mean that the model is an accurate 

predictor and the model may become over-fitted. Consequently, R~ dj is used as an indicator of model adequacy. 


K,j = i- 


SS E /(n — p) 
SS T /(n- 1) 


(18) 


The adjusted R~ penalizes the model for adding terms that are not meaningful, so it is most useful in evaluating and 
comparing candidate regression models, particularly when working with deterministic simulations. Cross validation 
is useful when collecting new data for validation purposes, but it is not possible. A reasonable procedure is to split 
the available data into two parts called the estimation data and the prediction data. 29 The estimation data are used to 
build the regression model, and the prediction data are then used to study the predictive ability of the model. 26 One 
useful statistic for cross validation is the prediction error sum of squares (PRESS). The PRESS residuals are defined 


as = y t — j where 5’(,) ' s the predicted value of the ; th observed response based on a model fit to the 
remaining n-1 sample points. 26 Thus, the PRESS statistic is defined as a measure of model quality; see Eq. (20). 

__ 2 z'' "\ 2 


PRESS = 


L y , - >%,_ 


i - /?., 


(19) 


Generally, when using the PRESS statistic, small values of PRESS are desired. Using the PRESS statistic, the 
predicted R 2 is described as 


r>2 

^ p red 


= 1 - 


PRESS 

SS-r 


( 20 ) 


This statistic gives an indication of the predictive capability of the regression model and a sense of how much 
variability in new observations the model might be expected to explain. 26 Table 3 summarizes the fit statistics. 


The final validation procedure is to apply a few confirmation runs using point prediction to test the 
regression metamodels. 26 A few points within the design space were selected that were not used to build the 
regression model. The results of the regression metamodel were compared with the simulation results. The percent 
difference was calculated to provide a measure of prediction for the regression metamodel, compared to the 
simulation. The percent difference between the metamodel and simulation ranged from 0.1% to 7%. 
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Table 3. Regression metamodel fit summary statistics. 



In-Phase Lift Coefficient 


Std. Dev. 

0.0235 

R-Squared 

0.8957 

Mean 

1.1047 

Adj R-Squared 

0.8305 



Pred R-Squared 

0.7670 


Out-of-phase Lift Coefficient 


Std. Dev. 

0.2294 

R-Squared 

0.9574 

Mean 

12.6622 

Adj R-Squared 

0.9307 



Pred R-Squared 

0.9110 

In-phase Pitching Moment Coefficient 

Std. Dev. 

0.00273 

R-Squared 

0.9452 

Mean 

0.20025 

Adj R-Squared 

0.9109 



Pred R-Squared 

0.8762 

Out-of-phase Pitching Moment Coefficient 

Std. Dev. 

0.0436 

R-Squared 

0.9575 

Mean 

1.5613 

Adj R-Squared 

0.9310 



Pred R-Squared 

0.9089 


V. Monte Carlo Simulations 

The study used two different Monte Carlo simulation approaches for the uncertainty analysis. The Appendix 
contains flowcharts that detail the two approaches. The two approaches are (1) a Taylor series based sensitivity and 
(2) a regression model based sensitivity. 31 

A. Taylor Series Based Sensitivity Analysis 

The first approach uses the classical Taylor series based sensitivity analysis in conjunction with a Monte Carlo 
simulation to determine the uncertainty. For a Taylor series based sensitivity approach, the partial derivatives of Eq. 
21 were computed from the regression (coded variables) models. The sensitivities are computed at a given nominal 
setting of the factors and correlations between factors are not considered, Eq. 22. 


1 

CD 

.- ol 

dC, 

8C, 

SJ EQ 

dC, 

oB eq 

dIS 

EQ 






O) 

£ 

1 


dIS J 


( 21 ) 


The algorithm sampled a pseudo-population with 100,000 trials. Uniform distributions were chosen for 5J E q and 
5B E q and normal distributions were chosen for 5BL and SIS. No distribution was used for 5k. Finally, the 
procedure computes the in-phase and out-of-phase lift and pitching moment coefficients and concludes with 
summary statistics and confidence intervals. 
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B. Direct Regression Model Approach 


( 22 ) 


(23) 


The second approach is a direct regression model-based Monte Carlo method. A pseudo-population of the factors 
(i.e. J eq , B eq , BL, and IS) was directly substituted into the regression models. A run of 100,000 Monte Carlo trials 
was performed. The samples for each factor are direct inputs for the metamodel regression models. A flowchart of 
the procedure is provided in the Appendix. Again, no distribution was used for the reduced frequency factor. 


C. Test Cases 


Several test cases were run and are summarized in Table 4 using both Monte Carlo simulation approaches 
described previously. The first test case represents the ideal factor settings for a test rig operating at a low reduced 
frequency. The ideal settings were based on reported information from NASA LaRC. The equivalent inertia and 
damping are set at their maximum limit. Backlash is set at its minimum and input saturation is not present. The 
reader should refer back to Table 1 for values of these settings. The second test case uses the factor settings at the 
center of the design space. The third test case uses factor settings at the low levels within the design space. It also 
assumes that input saturation is significant. The final two cases were chosen randomly within the design space. 
Confidence intervals were calculated using the method of reference 31 for skewed distributions.’ 1 


VI. Results and Discussion 

All of the test cases showed highly skewed histograms of the in-phase and out-of-phase lift and pitching moment 
coefficients. An example is shown in Figure 6 for the in-phase lift coefficient. This demonstrates the high degree of 
nonlinearity in the system. Tables (4) and (5) summarize the 95% probability confidence intervals for all test cases 
for each Monte Carlo simulation approach. 
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C La (bar) 

Figure 6. Histogram for in-phase lift coefficient. 


Table 4. Summary of 95% confidence intervals using Taylor series based sensitivity analysis. 


Test Case 

95% Confidence Interval 

Test Case 

95% Confidence Interval 

#1 

1.0816 <C L <1.1967 
14.2512 <C, <15.3689 
0.1942 < C ma < 0.2070 
1.8614 <C m < 2.0782 

#4 

1.0185 <C L <1.1910 
14.2323 <C, <15.4076 

^<7 

0.1853 <C„, <0.2062 
1.8731 <C,„ <2.1218 

q 

#2 

1.1225 <C L <1.1842 
12.7585 < C, <13.4310 

^<7 

0.2023 < C ma < 0.2096 
1.5772 <C m <1.7004 

q 

#5 

I. 1643 <C L <1.3211 

I I. 7270 <C, <12.7796 
0.2111 <C nia <0.2303 
1.3899 <C m <1.6111 

<7 

#3 

1.0714 <C L <1.2526 
12.7585 <C, <15.4326 

^<7 

0.1919 <C ma <0.2138 
1.8436 <C„, <2.1181 

q 
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Table 5. Summary of 95% confidence intervals using direct regression model method. 


Test Case 

95% Confidence Interval 

Test Case 

95% Confidence Interval 

#1 

0.7664 <C L <1.3756 
10.5639 <C, <16.9628 

^q 

0.2290 <C m <0.3538 

m a 

1.1743 <C m <2.3590 

#4 

0.7786 <C, <1.3573 

L -‘a 

10.4190 <C, <15.5613 
0.1568 < C m <0.2257 

rn a 

1.1576 <C,„ <2.1568 

q 

#2 

0.9644 <C L <1.1387 
12.3625 <C, <13.9079 

•7 

0.1832 <C„ <0.2042 

m a 

1.5027 <C m <1.8014 

#5 

0.8593 <C L <1.4391 

8.5762 < C, <14.1788 

% 

0.1734 <C,„ <0.2438 

rn a 

0.8004 <C m <1.8465 

rn q 

#3 

0.8238 <C L <1.5077 
8.7298 <C, <15.1489 

q 

0.1624 <C m <0.2452 

m a 

0.8564 <C m <2.0756 

m q 




The sensitivity analysis via a Taylor series (TS) expansion is simple, but its use has several drawbacks. It is 
based on computing derivatives and numerical differentiation can be difficult for complex functions. The TS 
expansion assumes the response is nearly linear over a small range. It will not always be easy to determine whether 
the approximations involved using these methods are acceptable. The regression model (RM) Monte Carlo method 
does not suffer from this problem because it can reach an arbitrary level of accuracy. In addition, the TS method is 
strictly based on the statistical moments of each parameter and does not directly incorporate the parameter 
probability distribution(s). 31 One of the major drawbacks is that the method lacks coupled variable effects. If an 
input variable is statistically correlated, one-factor-at-a-time sensitivity analysis does not address such correlations. 
Attempting to evaluate all the potential combinations of input parameters can be unmanageable. Also, high order 
estimates may be necessary to adequately address highly skewed distributions. 

The regression metamodel (RM) Monte Carlo method obtains high accuracy by using a sufficiently large number 
of runs. The RM Monte Carlo method can suffer from improper selection of probability distributions either from 
inadequate data or lack of understanding of the underlying physical process. 

The TS Monte Carlo method may be used to obtain preliminary answers. However, the RM Monte Carlo 
method is preferred when error propagation with a complex system is being studied because the method is easily 
implemented and generally applicable. 

The TS Monte Carlo method predicts that the response variability ranges from approximately 3% to 6% of the 
mean responses over all test cases. On the other hand, the RM model Monte Carlo method predicts that the response 
variability ranges from approximately 3% to 40% of the mean responses over all test cases. 

The most significant finding of the study is that researchers using different test rigs can have, even with the same 
aircraft model, significantly varied results due only to the differences in system dynamic parameters. Many 
researchers often assume that bias errors can simply be averaged out or that they are negligible. However, the 
parameters play a larger role in that the systematic errors interact with the aerodynamics, even with proper 
calibration of the forced oscillation test rig. The large variability predicted by the RM Monte Carlo method 
indicates that the bias effects are interacting with the aerodynamics. For example, input saturation causes the 
dynamic contributions of lift coefficient and pitching moment coefficient to have missing data at the peaks and 
troughs of the sine wave; an example is shown in Figure 7. This result has also been shown in experimental results; 
refer to the work of Kim et al.~ 2 The poor prediction of the stability derivatives is problematic because the analysis 
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used to compute the in-phase and out-of-phase coefficients assumes that the sine fidelity is high. The analysis used 
simply solves for the area under the curve to determine the in-phase and out-of-phase coefficients. Therefore, one 
cause of the poor prediction is due to poor sine fidelity caused by input saturation. 



Time (sec) 

Figure 7: Incomplete sinusoidal waveform of dynamic contribution of lift coefficient 

The bias effect of backlash is also highly problematic. Although some backlash is necessary in all geared 
systems to allow a certain amount of clearance between the components transmitting the motion under the load to 
avoid interference, wear, and excessive heat generation, it can be particularly challenging if the backlash is large 
(i.e. 8 arcmin or higher). Backlash is also not important in applications where there is no load reversal or the 
position after a reversal is not critical. However, forced oscillation wind tunnel testing requires precision 
positioning with frequent load reversal. Backlash can directly influence the positioning. When the gears are 
separated, the load is no longer under control, control is only retained on the motor side. When the gears again 
mate, backlash can introduce a nonlinear vibration effect that may be on the same order of magnitude as the 
aerodynamic effects being measured. Furthermore, the coupling effects between the systematic errors make it very 
difficult to decouple them from the aerodynamic measurement. Referring to Table 6, the coded (i.e. centered and 
scaled) regression metamodels can be used to examine the magnitude of the effects. Factor labels are those given in 
table 1. 

Reduced frequency is the most dominant contributor to the in-phase and out-of-phase components; it is shown in 
bold in Table 6. This result is required since the aerodynamic model developed was based on reduced frequency. 
Interestingly, the coupling effects are about the same order of magnitude. The coupling makes it difficult to 
determine which systematic error is dominating and should be eliminated. For example, looking at the in-phase 
pitching moment coefficient, the interaction between backlash and input saturation (i.e. D*E) has a similar 
magnitude effect compared to an interaction between equivalent damping and backlash (i.e. B*D). The backlash 
and input saturation interaction also have the same order effect as the interaction between equivalent damping and 
reduced frequency (i.e. B*C). No one effect is dominating. Overall combined contributions are affecting the final 
result making it difficult to decouple the systematic errors from the aerodynamics. 
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Table 6: Coded Regression Metamodels 


In-Phase Lift Coefficient 

Out-of-Phase Lift Coefficient 

In-Phase Pitching Moment 
Coefficient 

Out-of-Phase Pitching 
Moment Coefficient 

Cl 


c, 

A 


C m 


C "’.i 


1.108E+00 


1.264E+01 


2.005E-01 


1.555E+00 


2.525E-02 

* A 

-3.061 E-01 

* A 

3.524E-03 

* A 

-5.061 E-02 

* A 

2.443E-02 

* B 

-2.369E-01 

* B 

2.585E-03 

* B 

-4.668E-02 

* B 

1.140E-01 

*C 

-1.688E+00 

*C 

1.761 E-02 

*C 

-3.224E-01 

‘C 

4.892E-02 

* D 

-5.170E-01 

* D 

5.765E-03 

* D 

-9.549E-02 

* D 

2.172E-02 

* E 

-1 .208E-02 

* E 

3.197E-03 

* E 

-9.387E-03 

* E 

5.852E-03 

* A * B 

-7.487E-02 

* A* B 

7.959E-04 

* A * B 

-1 .309E-02 

* A * B 

-1.168E-02 

* A * E 

9.933E-02 

* A* E 

-1.210E-03 

* A * E 

2.023E-02 

* A* E 

1 .386E-02 

* B*C 

-1.373E-01 

* B * C 

1 .558E-03 

* B‘C 

-2.661 E-02 

* B * C 

-8.176E-03 

* B * D 

8.054E-02 

* B * D 

-9.214E-04 

* B * D 

1 .476E-02 

* B * D 

8.376E-03 

* C * D 

-7.797E-02 

* C * D 

8.947E-04 

* C * D 

-1.517E-02 

* C * D 

-8.135E-03 

* D * E 

9.299E-02 

* D* E 

-1.029E-03 

* D * E 

1 .690E-02 

* D* E 

1.215E-02 

* A A 2 

-1.1 12E-01 

* A A 2 

1 .544E-03 

* A A 2 

-2.066E-02 

* A A 2 

8.037E-03 

* B A 2 

-3.725E-02 

* B A 2 

7.871 E-04 

* B A 2 

-1.030E-02 

* B A 2 

4.221 E-02 

* C A 2 

-4.648E-01 

* C A 2 

5.321 E-03 

* C A 2 

-8.230E-02 

* C A 2 

-2.091 E-02 

* D A 2 

2.596E-01 

* D A 2 

-2.553E-03 

* D A 2 

4.557E-02 

* D A 2 

-5.023E-02 

* E A 2 

3.979E-01 

* E A 2 

-5.770E-03 

* E A 2 

8.170E-02 

* E A 2 

-3.757E-02 

* A A 3 

4.459E-01 

* A A 3 

-4.995E-03 

* A A 3 

7.620E-02 

* A A 3 

-2.540E-02 

* B A 3 

2.616E-01 

* B A 3 

-2.786E-03 

* B A 3 

5.075E-02 

* B A 3 

-5.476E-02 

* C A 3 

5.415E-01 

* C A 3 

-6.31 IE-03 

* C A 3 

1.038E-01 

* C A 3 

-4.631 E-02 

* D A 3 

4.968E-01 

* D A 3 

-5.529E-03 

* D A 3 

9.133E-02 

* D A 3 


VII. Conclusions 

Poor prediction of aircraft stability derivatives is often caused by unsteady, nonlinear aerodynamics brought on 
by high angle-of-attack and/or high angular rate maneuvers. Consequently, researchers have begun to use unsteady 
nonlinear modeling methodologies and forced oscillation wind tunnel testing in attempts to estimate dynamic 
stability derivatives. 

There have been measurement data inconsistencies noted by researchers. The inconsistencies are possibly from 
low accuracy and non-physical values of parameter estimation. Other sources of error are due to model structure 
error in the form of measurement bias. Bias errors are a particular concern because the aerodynamic phenomena are 
occurring below 10 Hz, which can interact with the test rig dynamics. Other sources of error not considered here 
include those inherent to all wind tunnel testing such as Reynolds number scaling, dynamic scaling, differences in 
turbulence intensity (transition), proximity to supports or tunnel boundaries. It should be noted that this study is 
only a foundation upon which to build. Many other sources of error may be considered, for example the force 
balance dynamics were not included in the study. 

This study showed variability in the responses in keeping with a system that is highly nonlinear. Responses were 
highly sensitive to the factors (i.e. J eq , B eq , k, IS, BL). In addition, and unfortunately, the nonlinear behavior of the 
dynamic system and its systematic biases has a coupling effect with the aerodynamics of an aircraft model (in this 
case the F-16XL). While factors may not change over the ranges chosen for this study in every facility, it is highly 
likely that similar factors will affect the results from forced oscillation testing. 
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Appendix 


Monte Carlo Simulation Based on Taylor Series Expansion: 



17 

American Institute of Aeronautics and Astronautics 








Monte Carlo Simulation Based on Regression Metamodels: 
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